      subroutine aysub(ay, aymin, aymax, bx, bz, iwid, jwid, kwid, ibp1, jbp1, &
         kbp1, dx, dy, dz) 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: iwid 
      integer , intent(in) :: jwid 
      integer , intent(in) :: kwid 
      integer , intent(in) :: ibp1 
      integer , intent(in) :: jbp1 
      integer , intent(in) :: kbp1 
      real , intent(out) :: aymin 
      real , intent(out) :: aymax 
      real , intent(in) :: dx 
      real  :: dy 
      real , intent(in) :: dz 
      real , intent(inout) :: ay(*) 
      real , intent(in) :: bx(*) 
      real , intent(in) :: bz(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: ibp2, jbp2, kbp2, j, ijk, i, k 
!-----------------------------------------------
      ibp2 = ibp1 + 1 
      jbp2 = jbp1 + 1 
      kbp2 = kbp1 + 1 
 
      do j = 1, jbp2 
         ijk = 1 + (2 - 1)*iwid + (j - 1)*jwid + (kbp2/2)*kwid 
         ay(ijk) = 0. 
         do i = 3, ibp2 
            ijk = 1 + (i - 1)*iwid + (j - 1)*jwid + (kbp2/2)*kwid 
            ay(ijk) = ay(ijk-iwid) + 0.5*(bz(ijk-iwid)+bz(ijk-iwid-kwid))*dx 
         end do 
         do i = 2, ibp2 
            if (i == 2) then 
               do k = 2 + kbp2/2, kbp2 
                  ijk = 1 + (i - 1)*iwid + (j - 1)*jwid + (k - 1)*kwid 
                  ay(ijk) = ay(ijk-kwid) - bx(ijk-kwid)*dz 
               end do 
            else 
               if (i == ibp2) then 
                  do k = 2 + kbp2/2, kbp2 
                     ijk = 1 + (i - 1)*iwid + (j - 1)*jwid + (k - 1)*kwid 
                     ay(ijk) = ay(ijk-kwid) - bx(ijk-iwid-kwid)*dz 
                  end do 
               else 
                  do k = 2 + kbp2/2, kbp2 
                     ijk = 1 + (i - 1)*iwid + (j - 1)*jwid + (k - 1)*kwid 
                     ay(ijk) = ay(ijk-kwid) - 0.5*(bx(ijk-kwid)+bx(ijk-iwid-&
                        kwid))*dz 
                  end do 
               endif 
            endif 
            if (i == 2) then 
               do k = kbp2/2, 2, -1 
                  ijk = 1 + (i - 1)*iwid + (j - 1)*jwid + (k - 1)*kwid 
                  ay(ijk) = ay(ijk+kwid) + bx(ijk)*dz 
               end do 
            else 
               if (i == ibp2) then 
                  do k = kbp2/2, 2, -1 
                     ijk = 1 + (i - 1)*iwid + (j - 1)*jwid + (k - 1)*kwid 
                     ay(ijk) = ay(ijk+kwid) + bx(ijk-iwid)*dz 
                  end do 
               else 
                  do k = kbp2/2, 2, -1 
                     ijk = 1 + (i - 1)*iwid + (j - 1)*jwid + (k - 1)*kwid 
                     ay(ijk) = ay(ijk+kwid) + 0.5*(bx(ijk)+bx(ijk-iwid))*dz 
                  end do 
               endif 
            endif 
         end do 
      end do 
      aymin = 1.D10 
      aymax = -1.D10 
      do i = 3, ibp2 
         ijk = 1 + (i - 1)*iwid + 1*jwid + (kbp2/2)*kwid 
         aymin = min(aymin,ay(ijk)) 
         aymax = max(aymax,ay(ijk)) 
      end do 
      return  
      end subroutine aysub 

! -------------------------------------------------------------------
!
      subroutine opoints(x,bzn,iwid,jwid,kwid,ibp1,jbp1,kbp1,rfsolar)
!
!	routine to compute distance between o points in 2 magnetic islands 
!	coalescence
!
      implicit none
      real(8) :: x(*),bzn(*),rfsolar
      integer :: ibp1,jbp1,kbp1,iwid,jwid,kwid,ibp2,jbp2,kbp2
!
      real(8) :: xext,xlast,xfirst
      integer :: i,j,k,ijk,imjk,istart,iend
! 
      ibp2 = ibp1 + 1 
      jbp2 = jbp1 + 1 
      kbp2 = kbp1 + 1 
      k=kbp2/2
      j=2
      istart=int(ibp2/5)
      iend=ibp2-istart
      
      xfirst=1d10
      xlast=-1d10
      do i=istart,iend
      ijk=1+(i-1)*iwid+(j-1)*jwid+(k-1)*kwid
      imjk=1+((i-1)-1)*iwid+(j-1)*jwid+(k-1)*kwid
      if(bzn(ijk)*bzn(imjk).le.0d0) then
!
!	x or o point found
!
      xext=x(imjk)-(x(ijk)-x(imjk))/(bzn(ijk)-bzn(imjk)+1d-10)*bzn(imjk)
      xfirst=min(xfirst,xext)
      xlast=max(xlast,xext)
      end if
      enddo
      rfsolar=xlast-xfirst
!      write(*,*)'rfsolar',rfsolar
      return
      end subroutine opoints
